vars =  {   'Coefftrue', 'Coefficients', 'StandardErrors', 'ComputationTime',...
    'OwnElast', 'CrossElastSameSeg', 'CrossElastDiffSeg',...
    'OneOwnElast','OneCrossElastSameSeg','OneCrossElastDiffSeg',...
    'DiversionSameSeg', 'DiversionDiffSeg', 'OneDiversionSameSeg','OneDiversionDiffSeg',...
    'BiasOwnElast', 'BiasCrossElastSameSeg', 'BiasCrossElastDiffSeg',...
    'FunctionValue', 'ExitFlag', 'AIC', 'BIC', ...
    'SaveSegment', 'SaveX1', 'OutsideShare', 'corrDstarX1'};

Nameresults = { '1resultsRCNL_LowCorr_Pho03Sigma1_Dstar0.mat';...
                '2resultsRCNL_LowCorr_Pho05Sigma05_Dstar0.mat';...
                '3resultsRCNL_HighCorr_Pho03Sigma1_Dstar0.mat';...
                '4resultsRCNL_HighCorr_Pho05Sigma05_Dstar0.mat';...
                '5resultsRCNL_LowCorr_Pho03Sigma1_Dstar1.mat';...
                '6resultsRCNL_LowCorr_Pho05Sigma05_Dstar1.mat';
                '7resultsRCNL_HighCorr_Pho03Sigma1_Dstar1.mat';...
                '8resultsRCNL_HighCorr_Pho05Sigma05_Dstar1.mat'};

TitleSheets = { '1LowCorr_Pho03Sgm1_Dstar0';...
                '2LowCorr_Pho05Sgm05_Dstar0';...
                '3HighCorr_Pho03Sgm1_Dstar0';...
                '4HighCorr_Pho05Sgm05_Dstar0';...
                '5LowCorr_Pho03Sgm1_Dstar1';...
                '6LowCorr_Pho05Sgm05_Dstar1';
                '7HighCorr_Pho03Sgm1_Dstar1';...
                '8HighCorr_Pho05Sgm05_Dstar1'};

File='RCNLresults';

for j = 1:8;
    
    load(Nameresults{j},vars{:});
    TitleSheet = TitleSheets{j};
    
    %% Adapt
    
    rho_names   = ['rc',num2str(Coefftrue(4,1))];
    rc_names    = ['sigma',num2str(Coefftrue(5,1))];
        
    %% Logit regression Dummy Segment and X1 (% correct classification)
    
    % define optimization options
    options = optimset('Display','off','FunValCheck','on','MaxFunEvals',Inf, ...
    'MaxIter',Inf,'TolFun',1e-5,'TolX',1e-6);
    ndraws      = 1000;
    pct         = zeros(ndraws,1);
    const       = ones(1250,1); % need a constant regressor too
    warning off
    
    for i = 1:ndraws
        % define logistic function
        logfun = @(x) 1./(1+exp(-x));
        % construct regressor matrix (data points x regressors)
        y  = SaveSegment(:,i);
        X1 = SaveX1(:,i);
        X  = [X1 const];
        % the cost function to minimize is the negative-log-likelihood for the binomial logit
        % here, the eps is a hack to ensure that the log returns a finite number
        costfun     = @(pp) -sum((y .* log(feval(logfun,X*pp')+eps)) + ...
            ((1-y) .* log(1-feval(logfun,X*pp')+eps)));
        x0          = zeros(1,size(X,2));
        params      = fminunc(costfun,x0,options);
        prediction  = X*params'>=0.5;
        % calculate percent correct
        pct(i) = sum(prediction==y) / length(y) * 100;
    end
    
    warning on
    %% Average classification
    AvCorrectClassified     = mean(pct);
    
    % Analysis of results
    
    %% Inspect results by visual look at the histogram of the coefficient distributions
    
    % Nested logit
%     figureName  = 'Distribution Pho';
%     figure2save = figure('Name',[num2str(rho_names),figureName],'Color',[1 1 1],'Position',get(0,'ScreenSize'));
    
%     histfit (squeeze((Coefficients(4,:,2,:))));
%     saveas(figure2save,[figureName,'.bmp']);
%     normplot(squeeze((Coefficients(4,:,2,:))));
    
    % Random Coefficient Logit
%     figureName  = 'Distribution Sigma';
%     figure2save = figure('Name',[num2str(rc_names),figureName],'Color',[1 1 1],'Position',get(0,'ScreenSize'));
    
%     histfit (squeeze((Coefficients(5,:,3,:))));
%     saveas(figure2save,[figureName,'.bmp']);
%     normplot(squeeze((Coefficients(5,:,3,:))));
    
    %% Coefficients and standard errors
    
    MeanCoefficient         = nanmean(Coefficients,4);
    EmpiricalStErrors       = nanstd(Coefficients,0,4);
    MeanStandardErrors      = nanmean(StandardErrors,4);
    MeanOutsideShare        = mean(OutsideShare,2);
    CompTime                = sum(ComputationTime,4);
    
    % Lilliefors Normality test
    % Nesting parameter
    [hNL,pNL]                 = lillietest(squeeze(Coefficients(4,:,2,:))); % Nested logit
    [hRCNLPho,pRCNLPho]       = lillietest(squeeze(Coefficients(4,:,4,:))); % RCNL
    % Random coefficient parameter
    [hRCL,pRCL]               = lillietest(squeeze(Coefficients(5,:,3,:))); % RCL
    [hRCNLSigma,pRCNLSigma]   = lillietest(squeeze(Coefficients(4,:,4,:))); % RCNL
    
    % Jarque-Bera Normality test
    % Nesting parameter
    [JBhNL,JBpNL]                 = jbtest(squeeze(Coefficients(4,:,2,:))); % Nested logit
    [JBhRCNLPho,JBpRCNLPho]       = jbtest(squeeze(Coefficients(4,:,4,:))); % RCNL
    % Random coefficient parameter
    [JBhRCL,JBpRCL]               = jbtest(squeeze(Coefficients(5,:,3,:))); % RCL
    [JBhRCNLSigma,JBpRCNLSigma]   = jbtest(squeeze(Coefficients(4,:,4,:))); % RCNL
    
    % Kolmogorov-Smirnov (KS) Test
    % Nesting parameter
    [KShNL,KSpNL]                 = kstest(squeeze(Coefficients(4,:,2,:))); % Nested logit
    [KShRCNLPho,KSpRCNLPho]       = kstest(squeeze(Coefficients(4,:,4,:))); % RCNL
    % Random coefficient parameter
    [KShRCL,KSpRCL]               = kstest(squeeze(Coefficients(5,:,3,:))); % RCL
    [KShRCNLSigma,KSpRCNLSigma]   = kstest(squeeze(Coefficients(4,:,4,:))); % RCNL
    
    
    %% AIC and BIC criteria to choose the appropriate model
    
    AICrc   = squeeze(AIC); % rows=number of estimators; cols=number of MC simulations
    MC      = size(AICrc,2);
    MinAIC  = zeros(1,MC);
    
    for j=1:MC
        [MinAIC(:,j)]     = find(AICrc(:,j)==min(AICrc(:,j)));
    end
    
    BICrc      = squeeze(BIC); % rows=number of estimators; cols=number of MC simulations
    MinBIC     = zeros(1,MC);
    
    for j=1:MC
        [MinBIC(:,j)]     = find(BICrc(:,j)==min(BICrc(:,j)));
    end
    
    CountAIC1 = countmember(1,MinAIC); % count how many times Logit model is the preferred model according to AIC criterion
    CountBIC1 = countmember(1,MinBIC); % count how many times Logit model is the preferred model according to BIC criterion
    
    CountAIC2 = countmember(2,MinAIC); % count how many times NL model is the preferred model according to AIC criterion
    CountBIC2 = countmember(2,MinBIC); % count how many times NL model is the preferred model according to BIC criterion
    
    CountAIC3 = countmember(3,MinAIC); % count how many times RCL model is the preferred model according to AIC criterion
    CountBIC3 = countmember(3,MinBIC); % count how many times RCL model is the preferred model according to BIC criterion
    
    CountAIC4 = countmember(4,MinAIC); % count how many times RCNL model is the preferred model according to AIC criterion
    CountBIC4 = countmember(4,MinBIC); % count how many times RCNL model is the preferred model according to BIC criterion
    
    %% Elasticity analysis
    
    MeanOwnElasticity         = nanmean(OwnElast,4);
    StDevOwnElasticity        = nanstd(OwnElast,0,4);
    
    MeanCrossElastSameSeg     = nanmean(CrossElastSameSeg,4);
    StDevCrossElastSameSeg    = nanstd(CrossElastSameSeg,0,4);
    
    MeanCrossElastDiffSeg     = nanmean(CrossElastDiffSeg,4);
    StDevCrossElastDiffSeg    = nanstd(CrossElastDiffSeg,0,4);

    OneOwnElasticity          = nanmean(OneOwnElast,4);
    OneStDevOwnElasticity     = nanstd(OneOwnElast,0,4);
    
    OneCrossElasticSameSeg    = nanmean(OneCrossElastSameSeg,4);
    StDevCrossElasticSameSeg  = nanstd(OneCrossElastSameSeg,0,4);
    
    OneCrossElasticDiffSeg    = nanmean(OneCrossElastDiffSeg,4);
    StDevCrossElasticDiffSeg  = nanstd(OneCrossElastDiffSeg,0,4);
   
    MeanDiversionSameSeg      = nanmean(DiversionSameSeg,4);
    StDevDiversionSameSeg     = nanstd(DiversionSameSeg,0,4);
    
    MeanDiversionDiffSeg      = nanmean(DiversionDiffSeg,4);
    StDevDiversionDiffSeg     = nanstd(DiversionDiffSeg,0,4);

    MeanOneDiversionSameSeg   = nanmean(OneDiversionSameSeg,4);
    StDevOneDiversionSameSeg  = nanstd(OneDiversionSameSeg,0,4);
        
    MeanOneDiversionDiffSeg   = nanmean(OneDiversionDiffSeg,4);
    StDevOneDiversionDiffSeg  = nanstd(OneDiversionDiffSeg,0,4);
    
    % Display and save results
    
    variables       = {'const';'dummy';'x1';'pho';'sigma'};
    titleCoeff      = {'Coefficients'};
    titleEstimator  = {'True','Logit','NL','RCL','RCNL'};
    resultsCoeff    = [Coefftrue, squeeze(MeanCoefficient)];
    
    disp(titleCoeff);
    disp([variables,num2cell(resultsCoeff)])
    
    variables       = {'const';'dummy';'x1';'pho';'sigma'};
    titleStErr1     = {'Empirical StError'};
    resultsStErr1   = squeeze(EmpiricalStErrors);
    
    disp(titleStErr1);
    disp([variables,num2cell(resultsStErr1)])
    
    variables       = {'const';'dummy';'x1';'pho';'sigma'};
    titleStErr2     = {'Calculated StError'};
    resultsStErr2   = squeeze(MeanStandardErrors);
    
    disp(titleStErr2);
    disp([variables,num2cell(resultsStErr2)])
    
    % Save results in excel
    
    xlswrite(File,titleCoeff, TitleSheet,'A1')
    xlswrite(File,titleEstimator,TitleSheet,'B1')
    xlswrite(File,resultsCoeff,TitleSheet,'B2')
    xlswrite(File,variables,TitleSheet,'A2')
    
    xlswrite(File,titleStErr1,TitleSheet,'A8')
    xlswrite(File,titleEstimator,TitleSheet,'B8')
    xlswrite(File,resultsStErr1,TitleSheet,'C9')
    xlswrite(File,variables,TitleSheet,'A9')
    
    xlswrite(File,titleStErr2,TitleSheet,'A15')
    xlswrite(File,titleEstimator,TitleSheet,'B15')
    xlswrite(File,resultsStErr2,TitleSheet,'C16')
    xlswrite(File,variables,TitleSheet,'A16')
    
    variables  = {'mins0';'avgs0';'maxs0'};
    xlswrite(File,variables,TitleSheet,'A23')
    xlswrite(File,squeeze(MeanOutsideShare),TitleSheet,'B23')
    
    variables  = {'CPUTime'};
    xlswrite(File,variables,TitleSheet,'A27')
    xlswrite(File,(squeeze(CompTime))',TitleSheet,'C27')
    
    % Elasticities
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'Own Elasticity'};
    
    xlswrite(File,titleElast,TitleSheet,'A29')
    xlswrite(File,variables,TitleSheet,'A30')
    xlswrite(File,squeeze(MeanOwnElasticity(:,1,:)),TitleSheet,'C30')
    
    titleElast = {'Cross Elasticities Same Seg'};
    
    xlswrite(File,titleElast,TitleSheet,'A33')
    xlswrite(File,variables,TitleSheet,'A34')
    xlswrite(File,squeeze(MeanCrossElastSameSeg(:,1,:)),TitleSheet,'C34')
    
    titleElast = {'Cross Elasticities Diff Seg'};
    
    xlswrite(File,titleElast,TitleSheet,'A37')
    xlswrite(File,variables,TitleSheet,'A38')
    xlswrite(File,squeeze(MeanCrossElastDiffSeg(:,1,:)),TitleSheet,'C38')
    
    %% AIC and BIC
    
    titleAIC    = {'AIC'};
    xlswrite(File,titleAIC,TitleSheet,'A41')
    xlswrite(File,titleEstimator,TitleSheet,'B41')
    xlswrite(File,[CountAIC1,CountAIC2,CountAIC3,CountAIC4],TitleSheet,'C42')
    
    titleBIC    = {'BIC'};
    xlswrite(File,titleBIC,TitleSheet,'A44')
    xlswrite(File,titleEstimator,TitleSheet,'B44')
    xlswrite(File,[CountBIC1,CountBIC2,CountBIC3,CountBIC4],TitleSheet,'C45')
    
    %% Lilliefors Normality tests
    
    title    = {'Lilliefors Normality test Nest'};
    xlswrite(File,title,TitleSheet,'A47')
    xlswrite(File,titleEstimator,TitleSheet,'B47')
    xlswrite(File,[NaN,hNL,NaN,hRCNLPho;NaN,pNL,NaN,pRCNLPho],TitleSheet,'C48')
    
    title    = {'Lilliefors Normality test RC'};
    xlswrite(File,title,TitleSheet,'A51')
    xlswrite(File,titleEstimator,TitleSheet,'B51')
    xlswrite(File,[NaN,NaN,hRCL,hRCNLSigma;NaN,NaN,pRCL,pRCNLSigma],TitleSheet,'C52')
    
    %% Jarque-Bera Normality tests
    
    title    = {'Jarque-Bera Normality test Nest'};
    xlswrite(File,title,TitleSheet,'A54')
    xlswrite(File,titleEstimator,TitleSheet,'B54')
    xlswrite(File,[NaN,hNL,NaN,JBhRCNLPho;NaN,pNL,NaN,JBpRCNLPho],TitleSheet,'C55')
    
    title    = {'Jarque-Bera Normality test RC'};
    xlswrite(File,title,TitleSheet,'A58')
    xlswrite(File,titleEstimator,TitleSheet,'B58')
    xlswrite(File,[NaN,NaN,JBhRCL,JBhRCNLSigma;NaN,NaN,JBpRCL,JBpRCNLSigma],TitleSheet,'C59')
    
    %% Kolmogorov-Smirnov Normality tests
    
    title    = {'Kolmogorov-Smirnov Normality test Nest'};
    xlswrite(File,title,TitleSheet,'A61')
    xlswrite(File,titleEstimator,TitleSheet,'B61')
    xlswrite(File,[NaN,hNL,NaN,KShRCNLPho;NaN,pNL,NaN,KSpRCNLPho],TitleSheet,'C62')
    
    title    = {'Kolmogorov-Smirnov Normality test RC'};
    xlswrite(File,title,TitleSheet,'A65')
    xlswrite(File,titleEstimator,TitleSheet,'B65')
    xlswrite(File,[NaN,NaN,KShRCL,KShRCNLSigma;NaN,NaN,KSpRCL,KSpRCNLSigma],TitleSheet,'C66')
    
    %% Diversion Ratios
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'Diversion Ratio Same Segment'};
    
    xlswrite(File,titleElast,TitleSheet,'A70')
    xlswrite(File,variables,TitleSheet,'A71')
    xlswrite(File,squeeze(MeanDiversionSameSeg(:,1,:)),TitleSheet,'C71')
    
    titleElast = {'Diversion Ratio Diff Segment'};
    xlswrite(File,titleElast,TitleSheet,'A74')
    xlswrite(File,variables,TitleSheet,'A75')
    xlswrite(File,squeeze(MeanDiversionDiffSeg(:,1,:)),TitleSheet,'C75')
    
    titleElast = {'One Diversion Ratio Same Segment'};
    xlswrite(File,titleElast,TitleSheet,'A78')
    xlswrite(File,(squeeze(MeanOneDiversionSameSeg(:,1,:)))',TitleSheet,'C79')
    titleElast = {'One Diversion Ratio Diff Segment'};
    xlswrite(File,titleElast,TitleSheet,'A81')
    xlswrite(File,(squeeze(MeanOneDiversionDiffSeg(:,1,:)))',TitleSheet,'C82')
    
    %% correctly classied
    xlswrite(File,{'correct classified'},TitleSheet,'A84')
    xlswrite(File,AvCorrectClassified,TitleSheet,'C84')
    
    %% Standard deviation elasticities
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'Standard deviation Own Elasticity'};
    
    xlswrite(File,titleElast,TitleSheet,'A87')
    xlswrite(File,variables,TitleSheet,'A88')
    xlswrite(File,squeeze(StDevOwnElasticity(:,1,:)),TitleSheet,'C88')
    
    titleElast = {'Standard deviation Cross Elasticities Same Seg'};
    
    xlswrite(File,titleElast,TitleSheet,'A91')
    xlswrite(File,variables,TitleSheet,'A92')
    xlswrite(File,squeeze(StDevCrossElastSameSeg(:,1,:)),TitleSheet,'C92')
    
    titleElast = {'Standard deviation Cross Elasticities Diff Seg'};
    
    xlswrite(File,titleElast,TitleSheet,'A94')
    xlswrite(File,variables,TitleSheet,'A95')
    xlswrite(File,squeeze(StDevCrossElastDiffSeg(:,1,:)),TitleSheet,'C95')
    
    %% Standard deviation diversion ratio
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'St.Dev. Diversion Same Segment'};
    
    xlswrite(File,titleElast,TitleSheet,'A98')
    xlswrite(File,variables,TitleSheet,'A99')
    xlswrite(File,squeeze(StDevDiversionSameSeg(:,1,:)),TitleSheet,'C99')
    
    titleElast = {'Standard deviation Diversion Ratio Diff Segment'};
    
    xlswrite(File,titleElast,TitleSheet,'A102')
    xlswrite(File,variables,TitleSheet,'A103')
    xlswrite(File,squeeze(StDevDiversionDiffSeg(:,1,:)),TitleSheet,'C103')
    
    titleElast = {'Standard deviation One Diversion Ratio Same Segment'};
    
    xlswrite(File,titleElast,TitleSheet,'A106')
    xlswrite(File,(squeeze(StDevOneDiversionSameSeg(:,1,:)))',TitleSheet,'C107')
    
    titleElast = {'Standard deviation One Diversion Ratio Diff Segment'};
    
    xlswrite(File,titleElast,TitleSheet,'A109')
    xlswrite(File,(squeeze(StDevOneDiversionDiffSeg(:,1,:)))',TitleSheet,'C110')   
    
    % One elasticities
   
    titleElast = {'One Own Elasticity'};
    xlswrite(File,titleElast,TitleSheet,'A113')
    xlswrite(File,(squeeze(OneOwnElasticity(:,1,:)))',TitleSheet,'C114')
    
    titleElast = {'One Cross Elasticity Same Segment'};
    xlswrite(File,titleElast,TitleSheet,'A116')
    xlswrite(File,(squeeze(OneCrossElasticSameSeg(:,1,:)))',TitleSheet,'C117')

    titleElast = {'One Cross Elasticity Diff Segment'};
    xlswrite(File,titleElast,TitleSheet,'A119')
    xlswrite(File,(squeeze(OneCrossElasticDiffSeg(:,1,:)))',TitleSheet,'C120')    
   
    % One standard deviations elasticities
    
    titleElast = {'Standard deviation One Own Elasticity'};
    
    xlswrite(File,titleElast,TitleSheet,'A123')
    xlswrite(File,(squeeze(OneStDevOwnElasticity(:,1,:)))',TitleSheet,'C124')
    
    titleElast = {'Standard deviation One Cross Elasticity Same Seg'};
    
    xlswrite(File,titleElast,TitleSheet,'A126')
    xlswrite(File,(squeeze(StDevCrossElasticSameSeg(:,1,:)))',TitleSheet,'C127')
    
    titleElast = {'Standard deviation One Cross Elasticity Diff Seg'};
    
    xlswrite(File,titleElast,TitleSheet,'A129')
    xlswrite(File,(squeeze(StDevCrossElasticDiffSeg(:,1,:)))',TitleSheet,'C130')  
    
end